library(geobr)
library(cowplot)
library(sf)
library(tidyverse)
library(viridis)
library(sfhotspot)
library(units)
library(lwgeom) # otimizar junção dos grids
data_PE <- read.csv("~/pear/data_PE.csv", dec = ",") # 262mil ocorrencias de todas as espécies via GBIF
adapta_pe<-read.csv("~/pear/DadosBiod.csv", sep = ";", dec = ",") # Dados do Adapta Brasil sobre Integridade dos Biomas
adapta_pe %>%
as_tibble() %>%
rename(code_muni=geocod_ibge) %>%
filter(code_muni!= "2605459") ->adapta_pePEAR - Biodiversidade & Ecossistemas
Respositório github para estas análises
Carregar dados e blibliotecas
Carregando bases espaciais dos municípios de PE
Usamos o pacote geobr para contruir os mapas que serão usados para ilustrar os resultados espaccialmente.
# # Correr os códigos abaixo somente uma vez
# muni_pe <- read_municipality(
# code_muni = "PE",
# year= 2022,
# showProgress = FALSE,
# simplified = FALSE
# )
# class(muni_pe)
#
# write_sf(muni_pe,"muni_pe.gpkg")
muni_pe<-st_read("muni_pe.gpkg")Reading layer `muni_pe' from data source `/home/felipe/pear/muni_pe.gpkg' using driver `GPKG'
Simple feature collection with 185 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -41.35834 ymin: -9.482897 xmax: -32.37777 ymax: -3.804622
Geodetic CRS: SIRGAS 2000
muni_pe<-muni_pe[muni_pe$code_muni != "2605459", ] # retirar Fernando de NoronhaGráfico teste
Um gráfico teste para unir os municípios de PE e os registros de espécies compilados no GBIF. Esses dados serão posteriormente usados para interpretar o risco à biodiversidade e ecossistemas e separar zonas mais ou menos amostradas quanto à biodiversidade.
# plot unindo registros do GBIF e
ggplot() +
geom_sf(data=muni_pe, fill="#2D3E50", color="#FEBF57", size=.15, show.legend = FALSE) +
coord_sf(xlim = c(-42, -34), ylim = c(-7, -10))+
geom_point(data = data_PE,aes(x=decimalLongitude, y=decimalLatitude),color = "red", alpha = 0.2)+
theme_minimal()## Registros por município
data_PE %>%
select(decimalLatitude,decimalLongitude) %>%
mutate(id_points = as.character(1:n()))->sp_coord
muni_pe <- st_as_sf(muni_pe, crs = 4326) # convert to sf object
sp_coord<-st_as_sf(data_PE, coords = c( "decimalLongitude", "decimalLatitude"), crs = 4326)
sp_coord <- st_transform(sp_coord, st_crs(muni_pe))
st_crs(sp_coord) # SIRGAS 2000Coordinate Reference System:
User input: SIRGAS 2000
wkt:
GEOGCRS["SIRGAS 2000",
DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["Latin America - Central America and South America - onshore and offshore. Brazil - onshore and offshore."],
BBOX[-59.87,-122.19,32.72,-25.28]],
ID["EPSG",4674]]
st_crs(muni_pe)# SIRGAS 2000 As duas OKCoordinate Reference System:
User input: SIRGAS 2000
wkt:
GEOGCRS["SIRGAS 2000",
DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["Latin America - Central America and South America - onshore and offshore. Brazil - onshore and offshore."],
BBOX[-59.87,-122.19,32.72,-25.28]],
ID["EPSG",4674]]
sp_coord$geometry # 262.332 pontos válidosGeometry set for 262332 features
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -41.33322 ymin: -9.474889 xmax: -32.38081 ymax: -3.807668
Geodetic CRS: SIRGAS 2000
First 5 geometries:
muni_pe$geom # 185 municípiosGeometry set for 184 features
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -41.35834 ymin: -9.482897 xmax: -34.80669 ymax: -7.272954
Geodetic CRS: SIRGAS 2000
First 5 geometries:
sp_counts_mun <- hotspot_count(sp_coord, grid = muni_pe) # 12,908 point is outside the area covered by the supplied polygons. Imprecisão dos datums?
sum(sp_counts_mun$n) # 258244 registros. É bastente[1] 249343
glimpse(sp_counts_mun) # OK para registrosRows: 184
Columns: 9
$ code_muni <dbl> 2600054, 2600104, 2600203, 2600302, 2600401, 2600500, 260…
$ name_muni <chr> "Abreu E Lima", "Afogados da Ingazeira", "Afrânio", "Agre…
$ code_state <dbl> 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 2…
$ abbrev_state <chr> "PE", "PE", "PE", "PE", "PE", "PE", "PE", "PE", "PE", "PE…
$ name_state <chr> "Pernambuco", "Pernambuco", "Pernambuco", "Pernambuco", "…
$ code_region <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ name_region <chr> "Nordeste", "Nordeste", "Nordeste", "Nordeste", "Nordeste…
$ n <int> 3371, 75, 843, 758, 734, 398, 147, 5286, 560, 135, 7, 247…
$ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-34.91405 -..., MULTIPOLYGON…
sp_counts_mun %>%
mutate(area_muni_km2=st_area(sp_counts_mun$geometry)/1000000) %>% # calcular área dos municípios
mutate(sp_per_area=n/area_muni_km2) %>% # espécies por KM²
filter(code_muni!= "2605459")->sp_counts_mun
sp_counts_mun %>%
mutate(rds=adapta_pe$rds) %>%
mutate(sp_per_area=as.numeric(sp_per_area))->sp_counts_munEsforço de coleta por RDS de PE
Temos, no geral um esforço de coleta baixo para PE em geral, com poucos registros por município e a região metropolitana (MET) é a mais bem amostrada. Destaque para o SFR, pouquíssimo amostrada na média.
sp_counts_mun %>%
as_tibble() %>%
left_join(adapta_pe[,-3], by="code_muni") %>%
#group_by(rds) %>%
#summarise(media = mean(sp_per_area), sd = sd(sp_per_area)) %>%
#arrange(desc(media)) %>%
#mutate(rds = fct_reorder(rds, mean(sp_per_area))) %>%
ggplot(aes(rds,log(sp_per_area)))+
geom_jitter(size=1)+
geom_errorbar(stat = "summary", fun.data = mean_cl_boot, width = 0.4)+
geom_point(size = 4, stat = "summary", fun = mean, shape = 1) +
#geom_errorbar(aes(ymin=media-sd, ymax=media+sd), width=.5,
# position=position_dodge(0.05))+
labs(x= "Região de Desenvolvimento Socioeconômico", y="Número médio de registros de espécies / km²")+
geom_hline(yintercept = mean(log(sp_counts_mun$sp_per_area)), col = "red", linetype="dashed")+
theme_bw()Dados de registros de espécies por município (fonte: GBIF)
# Gráficos e Mapas com a base completa biodiversidade
ggplot()+
geom_sf(data=sp_counts_mun, aes(fill=log(n+1), show.legend = TRUE)) +
scale_fill_viridis(name="")+
labs(title = "Número de registros de espécies por km²")Análise de risco climático para área temática de Biodievrsidade e Ecossistemas
Dados provenientes do Adapta Brasil
dados_comp_geom<-left_join(muni_pe[,c(1,8)],adapta_pe, by="code_muni")
dados_comp_geom %>%
st_transform(., crs=4674)Simple feature collection with 184 features and 36 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -41.35834 ymin: -9.482897 xmax: -34.80669 ymax: -7.272954
Geodetic CRS: SIRGAS 2000
First 10 features:
code_muni id rds nome vuln sensib s_luc_2017 sl_agrotox
1 2600054 4202 MET Abreu e Lima 0.44 0.33 0.34 0.00
2 2600104 4203 PAJ Afogados da Ingazeira 0.68 0.36 0.15 0.04
3 2600203 4204 SFR Afrânio 0.76 0.50 0.84 0.89
4 2600302 4205 AGC Agrestina 0.72 0.48 0.20 0.10
5 2600401 4206 MLS Água Preta 0.45 0.39 0.27 0.22
6 2600500 4207 AGM Águas Belas 0.68 0.40 0.14 0.05
7 2600609 4208 AGC Alagoinha 0.67 0.34 0.15 0.01
8 2600708 4209 MNS Aliança 0.64 0.35 0.32 0.08
9 2600807 4210 AGC Altinho 0.73 0.41 0.17 0.08
10 2600906 4211 MLS Amaraji 0.68 0.34 0.24 0.27
sl_agr_fam sl_prat_sust sl_fogo s_csol_2017 sc_passivo sc_desmat sc_pastag
1 0.68 0.96 0.00 0.07 0.00 0.19 0.13
2 0.15 0.75 0.14 0.18 0.06 0.30 0.83
3 0.45 0.84 0.85 0.24 0.10 0.39 0.98
4 0.30 0.71 0.03 0.28 0.17 0.59 0.39
5 0.32 0.79 0.10 0.08 0.76 0.08 0.08
6 0.23 0.67 0.05 0.28 0.17 0.38 0.91
7 0.29 0.84 0.02 0.20 0.02 0.38 0.95
8 0.20 0.54 0.44 0.17 0.96 0.04 0.29
9 0.25 0.84 0.07 0.32 0.12 0.52 0.59
10 0.20 0.93 0.07 0.12 0.22 0.16 0.09
sc_miner s_dens_pop_faz sd_popul sd_fazend s_infra_2017 si_rodov si_hidro
1 0 0.99 1.00 0.98 0.33 0.39 0.27
2 0 0.76 0.80 0.73 0.55 0.45 0.64
3 0 0.34 0.11 0.58 0.65 0.50 0.81
4 0 0.85 1.00 0.71 0.67 0.82 0.52
5 0 0.61 0.57 0.65 0.51 0.57 0.45
6 0 0.53 0.40 0.66 0.57 0.34 0.79
7 0 0.65 0.55 0.75 0.49 0.31 0.66
8 0 0.81 1.00 0.61 0.38 0.51 0.26
9 0 0.58 0.41 0.75 0.49 0.42 0.55
10 0 0.76 0.79 0.73 0.44 0.46 0.41
s_cap sc_orient sc_ucons expos e_cobexpo e_cobnat e_areaprot e_priocons
1 0.42 0.01 0.83 0.48 0.39 0.33 0.53 0.39
2 0.05 0.05 NA 0.53 0.49 0.49 1.00 0.32
3 0.04 0.04 NA 0.49 0.42 0.37 1.00 0.42
4 0.08 0.08 NA 0.74 0.88 0.88 1.00 0.42
5 0.46 0.46 NA 0.71 0.82 0.82 1.00 0.35
6 0.08 0.08 NA 0.62 0.65 0.68 0.00 0.65
7 0.05 0.05 NA 0.39 0.23 0.50 0.23 0.00
8 0.09 0.09 NA 0.80 0.99 0.97 1.00 0.99
9 0.00 0.01 0.00 0.68 0.78 0.78 1.00 0.23
10 0.02 0.02 NA 0.80 1.00 0.83 1.00 1.00
ameaca_atual ameaca_..SWL1.5 ameaca_SWL2.0 risco_atual risco_SWL1.5
1 0.39 0.04 0.25 0.58 0.10
2 0.01 0.24 0.24 0.03 0.59
3 0.02 0.24 0.24 0.05 0.60
4 0.41 0.45 0.53 0.78 0.80
5 0.47 0.08 0.24 0.70 0.35
6 0.09 0.14 0.13 0.44 0.52
7 0.03 0.25 0.26 0.09 0.54
8 0.67 0.23 0.24 0.87 0.66
9 0.37 0.45 0.47 0.75 0.79
10 0.51 0.20 0.18 0.83 0.64
risco_SWL2.0 geom
1 0.49 MULTIPOLYGON (((-34.91405 -...
2 0.59 MULTIPOLYGON (((-37.61725 -...
3 0.61 MULTIPOLYGON (((-41.03812 -...
4 0.83 MULTIPOLYGON (((-35.93697 -...
5 0.57 MULTIPOLYGON (((-35.38163 -...
6 0.50 MULTIPOLYGON (((-37.09861 -...
7 0.54 MULTIPOLYGON (((-36.78101 -...
8 0.67 MULTIPOLYGON (((-35.20328 -...
9 0.80 MULTIPOLYGON (((-36.09605 -...
10 0.62 MULTIPOLYGON (((-35.46141 -...
dados_comp_geom %>%
group_by(rds) %>%
summarise(across(geom, ~ st_combine(.)), .groups = "keep")->geom_rds # Funde bases de forma que as categorias de RDS são incluídas numa única base
# ggplot(data=geom_rds)+geom_sf(aes(fill=rds)) # só checando se está correto. OK
centroids_rds <- geom_rds %>%
st_make_valid() %>%
st_centroid() Análise de Vulnerabilidade
Grau de suscetibilidade do bioma a uma determinada ameaça climática. A vulnerabilidade refere-se às características que um determinado ecossistema possui que refletem à suscetibilidade dele em relação às variações climáticas que podem afetar, direta ou indiretamente, a integridade do bioma. A vulnerabilidade está vinculada a situação de sensibilidade e capacidade adaptativa desse sistema, ou seja, se baseia em fatores conjunturais existentes antes do impacto climático. (fonte: AdaptaBrasil)
Sensibilidade do bioma
Grau de alteração potencial do estado que o ecossistema pode sofrer, direta ou indiretamente, uma vez em contato com as variações climáticas que afetam a integridade do bioma. A sensibilidade é uma propriedade inerente a um ecossistema, existente antes da ameaça climática e independente (separada) da exposição. O Índice de Sensibilidade é composto pelos indicadores temáticos: Mudança e uso do solo, Cobertura do solo, Densidade populacional e de estabelecimentos, e Infraestrutura.
## Sensibilidade
library(RColorBrewer)
dados_comp_geom %>%
mutate(categ=case_when(
sensib <0.2 ~ "Muito Baixo",
sensib >=0.2 & sensib < 0.4 ~ "Baixo",
sensib >=0.4 & sensib < 0.6 ~ "Médio",
sensib >=0.6 & sensib < 0.8 ~ "Alto",
sensib >=0.8 ~ "Muito Alto",
TRUE ~ "unknown")) %>%
ggplot()+geom_sf(aes(fill=categ), show.legend = TRUE)+
scale_fill_brewer(palette="RdYlGn", name = "Categoria")+
labs(title = "Sensibilidade do Bioma")+
theme_bw()## PCA para entender componetes da sensibiidade por categoria de risco
dados_comp_geom %>%
as_tibble() %>%
mutate(categ=case_when(
sensib <0.2 ~ "Muito Baixo",
sensib >=0.2 & sensib < 0.4 ~ "Baixo",
sensib >=0.4 & sensib < 0.6 ~ "Médio",
sensib >=0.6 & sensib < 0.8 ~ "Alto",
sensib >=0.8 ~ "Muito Alto",
TRUE ~ "unknown")) %>%
select(starts_with("s_")|categ |rds)-> dados_sensib
pca_sensib<-prcomp(dados_sensib[,-c(5,6,7)], scale = TRUE)
summary(pca_sensib)Importance of components:
PC1 PC2 PC3 PC4
Standard deviation 1.2271 0.9893 0.9737 0.7533
Proportion of Variance 0.3764 0.2447 0.2370 0.1419
Cumulative Proportion 0.3764 0.6211 0.8581 1.0000
pca_sensib$rotation PC1 PC2 PC3 PC4
s_luc_2017 -0.5622488 -0.29858856 -0.5107138 0.5778343
s_csol_2017 -0.3018989 -0.51688230 0.7892945 0.1367627
s_dens_pop_faz 0.6751392 -0.04063902 0.1053063 0.7290035
s_infra_2017 -0.3700276 0.80126527 0.3241855 0.3405251
pca_data_sensib<-as.data.frame(pca_sensib$x[, 1:2]) # extract first two PCs
pca_data_sensib <- cbind(pca_data_sensib, dados_sensib$categ, dados_sensib$rds) # add species to df
colnames(pca_data_sensib) <- c("PC1", "PC2", "categoria","RDS") # change column names
pca_data_sensib %>%
ggplot() +
aes(PC1, PC2, color = categoria, shape = categoria) + # define plot area
geom_point(size = 2) + # adding data points
coord_fixed()+ # fixing coordinates
stat_ellipse(geom="polygon", level=0.95, alpha=0.1) #adding a statResumo da análise de sensibilidade
Os componentes da sensibilidade são praticamente ortogonais, significando que qualquer redução na dimensionalidade é prejudicial à interpretação. O primeiro PC é o úico com valor > 1 mas os PCs 2 e 3 tem valores muito próximos de 1. Portanto, o melhor é proceder para uma análise dos componentes do índice de sensibilidade de forma separada e isolada.
Sensibilidade do bioma
Daqui por diante as análises virão com duas representações gráficoas: 1) mapa de de PE com valores por município, podendo ser valores brutos ou categóricos (categorias Adapta Brasil) e 2) Gráfico mostrando valores médios e variação por região de desenvolvimento socioeconômico (RDS).
## RDS Mapa
col.map<-c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a','#ffff99','#b15928')
geom_rds %>%
ggplot()+geom_sf(aes(fill=rds))+
scale_fill_manual(values= col.map)+
geom_sf_label(aes(label = rds))+
theme_bw()+
theme(legend.position = "none")+
xlab("")+ylab("")->rds_map
rds_mapEssas são as RDS de Pernambuco
Indicador temático - Mudança de uso do solo
dados_comp_geom %>%
mutate(categ=case_when(
s_luc_2017 <0.2 ~ "Muito Baixo",
s_luc_2017 >=0.2 & s_luc_2017 < 0.4 ~ "Baixo",
s_luc_2017 >=0.4 & s_luc_2017 < 0.6 ~ "Médio",
s_luc_2017 >=0.6 & s_luc_2017 < 0.8 ~ "Alto",
s_luc_2017 >=0.8 ~ "Muito Alto")) %>%
arrange(desc(s_luc_2017)) %>%
mutate(categ = fct_reorder(categ, desc(s_luc_2017))) %>%
ggplot()+geom_sf(aes(fill=categ), show.legend = TRUE)+
scale_fill_brewer(palette="RdYlGn", name = "")+theme_bw()+theme(legend.position = "bottom")+
labs(title = "Sensibilidade à Mudanças no Uso de Solo")->luc_map
luc_map## Análise por RDS
dados_comp_geom %>%
mutate(categ=case_when(
s_luc_2017 <0.2 ~ "Muito Baixo",
s_luc_2017 >=0.2 & s_luc_2017 < 0.4 ~ "Baixo",
s_luc_2017 >=0.4 & s_luc_2017 < 0.6 ~ "Médio",
s_luc_2017 >=0.6 & s_luc_2017 < 0.8 ~ "Alto",
s_luc_2017 >=0.8 ~ "Muito Alto")) %>%
arrange(desc(s_luc_2017)) %>%
mutate(categ = fct_reorder(categ, desc(s_luc_2017))) %>%
ggplot(aes(rds,s_luc_2017))+
geom_jitter(size=1)+
geom_errorbar(stat = "summary", fun.data = mean_cl_boot, width = 0.4)+
geom_point(size = 4, stat = "summary", fun = mean, shape = 1) +
#geom_errorbar(aes(ymin=media-sd, ymax=media+sd), width=.5,
# position=position_dodge(0.05))+
labs(x= "Região de Desenvolvimento \n Socioeconômico", y="Sensibilidade à Mudanças no Uso de Solo")+
geom_hline(yintercept = mean(dados_comp_geom$s_luc_2017), col = "red", linetype="dashed")+
theme_bw()+coord_flip()->luc_graph
luc_graphplot_grid(luc_map, luc_graph, ncol=1, align="v", rel_widths = c(1,0.6), labels = c("a","b"))Valores médios de Vulnerabilidae pro RDS
\[ Vulnerabilidade = Sensibilidade+Exposição-Capacidade adaptativa \]
Grau de suscetibilidade do bioma a uma determinada ameaça climática. A vulnerabilidade refere-se às características que um determinado ecossistema possui que refletem à suscetibilidade dele em relação às variações climáticas que podem afetar, direta ou indiretamente, a integridade do bioma. A vulnerabilidade está vinculada a situação de sensibilidade e capacidade adaptativa desse sistema, ou seja, se baseia em fatores conjunturais existentes antes do impacto climático. (Fonte: Adapta Brasil)
## Todos os sub-índices
dados_comp_geom %>%
as_tibble() %>%
select(rds, sensib, s_cap, expos ) %>%
pivot_longer(!rds, names_to = "indic") %>%
group_by(rds, indic) %>%
summarise_if(is.numeric, list(mean, sd)) %>%
ggplot(aes(rds,fn1, color=indic))+
geom_point(size=3,position=position_dodge(width=0.7))+
scale_color_discrete(name= "",labels=c("Exposição", "Sensibilidade", "Capacidade Adaptativa"))+
geom_pointrange(aes(ymin=fn1-fn2, ymax=fn1+fn2),position=position_dodge(width=0.7))+
coord_flip()+
theme_bw()+
theme(legend.position = "top")+
labs(y="Valor do Índice", x="Região de Desenvolvimento Socioeconômico")->A
ANota-se aqui vários padrões:
Exposição é bastante variável entre e dentro das RDS, sugerindo que os municípios possuem contextos muito distintos que aumentam ou diminuem sua exposição ao risco climático
Sensibidade possui menor variação entre RDS mas algumas delas são compostas por municípios muito discrepantes (p.ex. MET e MLS )
A capacidade adapttiva é bastente uniforme entre as RDS com menor variação entre e dentro dos grupos, provavelmente indicando poucas variações nos indicadores individuais.
Risco Climático
## Risco climático atual e projetado
dados_comp_geom %>%
as_tibble() %>%
select(rds,risco_atual:risco_SWL2.0) %>%
pivot_longer(!rds, names_to = "risco") %>%
group_by(rds, risco) %>%
summarise_if(is.numeric, list(mean, sd)) %>%
ggplot(aes(rds,fn1, group=risco, color=risco))+
geom_point(size=3,position=position_dodge(width=0.7))+
scale_color_discrete(name= "",labels=c("Risco Atual", "Risco SWL 1.5", "Risco SWL 2.0"))+
geom_pointrange(aes(ymin=fn1-fn2, ymax=fn1+fn2),position=position_dodge(width=0.7))+
coord_flip()+
theme(legend.position = "top")+
labs(y="Risco Climático", x="Região de Desenvolvimento Socioeconômico")->B
B## Ameaça à Integridade do Bioma
dados_comp_geom %>%
as_tibble() %>%
select(rds,ameaca_atual:ameaca_SWL2.0) %>%
pivot_longer(!rds, names_to = "risco") %>%
group_by(rds, risco) %>%
summarise_if(is.numeric, list(mean, sd)) %>%
ggplot(aes(rds,fn1, group=risco, color=risco))+
geom_point(size=3,position=position_dodge(width=0.7))+
scale_color_discrete(name= "",labels=c("Ameaça Atual", "Ameaça SWL 1.5", "Ameaça SWL 2.0"))+
geom_pointrange(aes(ymin=fn1-fn2, ymax=fn1+fn2),position=position_dodge(width=0.7))+
coord_flip()+
theme(legend.position = "top")+
labs(y="Ameaça à Integridade do Bioma", x="Região de Desenvolvimento Socioeconômico")->C
C